#### This script produces the table with the models from the original analyses
#### using a different specification of the dependent variable (Gender Turnout Gap).
#### All code lines are directly copied and pasted from the original dataverse files 
#### by Skorge (2021) and modified to include the new dependent variable. 



##########################################################################
###            This part of the code replicates the models 4           ###
###           and 6 in Table 1 (p. 9) presented by the author in the   ###
###           paper's main analysis, but using Gender Turnout          ###
###                     Gap as the dependent variable.                 ###
##########################################################################

dv2 <- "turnoutGap"

# PR, FEs and all covariates
m.2a <- lm_robust(as.formula(paste0(dv2," ~ prFrom10 +",
                                    paste0(c(v.ctrls.pol,
                                             v.ctrls.other),collapse="+"))),
                  data=dt.a,
                  clusters=muni,
                  fixed_effects=~muni + year,
                  se_type = "stata")




# PR, FEs, all covariates, and municipality-specific time trends squared
m.2b <- lm_robust(as.formula(paste0(dv2," ~ prFrom10 +",
                                    paste0(c(v.ctrls.pol,
                                             v.ctrls.other,
                                             v_time,
                                             v_time_sq),collapse="+"))),
                  data=dt.a,
                  clusters=muni,
                  fixed_effects=~muni + year,
                  se_type = "stata")



##########################################################################
###  This part of the code replicates the model 6 in Table E.1 and     ###
###  model 4 in Table G.1 from the original Supplementary Information  ###
###                     provided by the author                         ###
##########################################################################

#####################################
### Model 6 Table E.1 (p. 7 - SI) ###
#####################################



# Using only data for 1916 to 1919 and municipalities present in 1916 with controls
m.2c <- lm_robust(as.formula(paste0(dv2," ~ prFrom10n +",
                                    paste0(v.ctrls,collapse="+"))),
                  data=dt,
                  subset=presentIn1916==1 & year %in% 1916:1919,
                  clusters=muni,
                  fixed_effects=~muni + year,
                  se_type="stata")

v.mean2.nc <- dt %>% 
  slice(as.numeric(attr(m.2c[["fitted.values"]],"names"))) %>% 
  filter(year>=1919 & pr16==0) %>%
  summarize(Mean(turnoutGap)) %>% 
  pull()



#####################################
### Model 4 Table G.1 (p. 7 - SI) ###
#####################################


# Model with leads and all covariates
m.2d <- lm_robust(as.formula(paste0("turnoutGap ~ prFrom10 +",
                                    paste0(c(v.leads,
                                             v.ctrls.pol,
                                             v.ctrls.other),collapse="+"))),
                  data=dt.a,
                  clusters=muni,
                  fixed_effects=~muni + year,
                  se_type="stata")

v.mean2.nd <- dt %>% 
  slice(as.numeric(attr(m.1d[["fitted.values"]],"names"))) %>% 
  filter(year>=1919 & pr16==0) %>%
  summarize(Mean(turnoutGap)) %>% 
  pull()

#### This chunk creates Table with the models replicated above

# Put all the models in the list
ls.m2 <- mget(ls(pattern="^m\\.2"))



# Use the get.regtable function to create
# regression output table


dt.reg2 <- get.regtable(ls.m2,
                       var.names = c("PR reform (1919)", "PR (lead t+1)","PR (lead t+2)"),
                       keep.vars = c("pr", "f1.prFrom10", "f2.prFrom10"),
                       digits.b.0=3,
                       weights=FALSE,
                       checkmark=c("Yes",""),
                       dv.mean=as.list(c(rep(dt.a[, Mean(turnoutGap[year>=1919 & pr16==0])],2),v.mean2.nc, v.mean2.nd)),
                       dv.mean.name="Mean $Y_{treated, t \\geq 1919}$",
                       effects="se",
                       model.labs=rep(c("Model 4, Table 1","Model 6, Table 1","Model 6, Table E.1","Model 4, Table G.1"),1))

# Add description of covariates included
dt.reg2[dt.reg$term %in% "Covariates",] <- c("Covariates", "All covars./No time trends","All covars. and time trends","All covars.","All covars. and Leads")


options(knitr.kable.NA = '')

kable.out <- dt.reg2[-c(1,12:14),] %>%
  kable(caption="OLS regression results reporting the estimated effect of PR on gender equality in voting, 1910-1928, with alternative DV (Gender Turnout Gap)",
        booktabs=TRUE,longtable=FALSE,format="html",
        align=c("l",rep("c",ncol(dt.reg)-1)),
        escape=FALSE,
        linesep="",
        row.names=FALSE,
        col.names=dt.reg[1,]) %>%
  kable_styling(font_size = 10,
                position = "center",
                full_width = FALSE) %>%
  row_spec(2, extra_latex_after="\\midrule") %>%
  row_spec(5, extra_latex_after="\\midrule") %>%
  add_header_above(c(" " = 1,
                     "Gender Turnout Gap" = 4
  )) %>%
  footnote(general='Standard errors clustered by municipality are in parentheses. These models are estimated identically to the author´s original models but using Gender Turnout Gap as the dependent variable. Gender Turnout Gap is mesured as measured as the difference between women’s and men’s turnout, i.e., [((Women votes cast)/(Eligible women to vote))×100]-[((Men votes cast)/(Eligible men to vote))×100]. Model 4, Table 1 reports the model in the main analysis including all the demographic and political covariates but no time trends; Model 6, Table 1 reports the same model including the linear and squared time trends (time trends are municipality specific); Model 6, Table E.1, reports the robustness check presented by the author in the SI where the full sample of  municipalities is used but time is restricted to the 1916 and 1919 elections (as in a more traditional Diff-in-Diff setting) using the full set of demographic and political covariates; Model 4, Table G.1 reports the Model 4 from the main analysis (All covars./No time trends) with leads of the 1919 PR treatment. The demographic and political covariates  included are: population (log), the female fraction of the population, eligible voters (log), the female fraction of eligible voters, and representatives in the municipal council (log), percentage of the population in nonconformist (Dissenting) religious societies and variables for the percentage of the working population in each of the employment categories industry, shipping, services, and agriculture (with agriculture as the omitted category). The effect size is the percentage increase from the counterfactual outcome: $\\\\hat{\\\\gamma}_{pr}/(\\\\bar{Y}_{treated, t \\\\geq 1919} - \\\\hat{\\\\gamma}_{pr}) \\\\times 100$, where $\\\\hat{\\\\gamma}_{pr}$ is the estimated PR coefficient and  $\\\\bar{Y}_{treated, t \\\\geq 1919}$ is the mean outcome in year $\\\\geq$ 1919 for the treated municipalities.',
           general_title="Note:",
           footnote_as_chunk=TRUE,
           threeparttable=TRUE,
           escape=FALSE)

kable.out

save_kable(kable.out, file=here("tables","table_replication_alternativeDV.html"))

setDT(dt.reg2)
